/*
 * @Description: 
 * @Autor: kaikaima
 * @Date: 2021-05-17 13:53:42
 * @LastEditors: kaikaima
 * @LastEditTime: 2021-05-25 21:02:31
 */

#include<iostream>
#include<cmath>

#ifndef __hhxxttxs_ITERATION__
#define __hhxxttxs_ITERATION__

class G_S {
    protected:
    unsigned N;
    double* A;
    double* x;
    double* B;
    double Up_Limit_Of_Error;
    double residual;
    public:
    G_S(unsigned _n) {
        N=_n;
        A=new double[N*N];
        x=new double[N];
        B=new double[N];
        Up_Limit_Of_Error=__DBL_MAX__;//avoid uninitialized data
        residual=0;
    };
    ~G_S() {
        delete[] A;
        delete[] x;
        delete[] B;
    };
    void set(double **_a){
        for(int i=0;i<N;i++)
        for(int j=0;j<N;j++)
            A[i*N+j]=_a[i][j];
        for(int i=0;i<N;i++)
        for(int j=0;j<i;j++){
            A[i*N+j]/=(-A[i*N+i]);
            A[j*N+i]/=(-A[j*N+j]);
        }
    };
    void set(double* _a){
        for(int i=0;i<N;i++){
            x[i]=B[i]=_a[i]/A[i*N+i];
        }
        residual=0;
        for(int i=0;i<N;i++){
            double key=x[i];
            x[i]=B[i];
            for(int j=0;j<N;j++)
                if(i!=j) x[i]+=(A[i*N+j]*x[j]);
            if(residual<fabs(x[i]-key)) residual=fabs(x[i]-key);
        }
    };
    void set(double _a){
        if(_a>0) Up_Limit_Of_Error=_a;
    };
    double* iteration(){
        short flag=1; //avoid endless loop
        while(residual>=Up_Limit_Of_Error && flag>0){
            flag++;
            residual=0;
            for(int i=0;i<N;i++){
                double key=x[i];
                x[i]=B[i];
                for(int j=0;j<N;j++)
                    if(i!=j) x[i]+=(A[i*N+j]*x[j]);
                if(residual<fabs(x[i]-key)) residual=fabs(x[i]-key);
            }
        }
        std::cout<<"flag="<<flag<<std::endl;
        return x;
    };
};

class Jacobi : public G_S {
    protected:
    double* y;
    public:
    Jacobi(unsigned _n) : G_S(_n) {
        y=new double [_n];
    };
    ~Jacobi(){
        delete[] y;
    };
    void set(double** _a){
        G_S::set(_a);
    };
    void set(double* _a){
        for(int i=0;i<N;i++){
            x[i]=B[i]=_a[i]/A[i*N+i];
        }
        residual=0;
        for(int i=0;i<N;i++){
            y[i]=B[i];
            for(int j=0;j<N;j++){
                if(i!=j) y[i]+=(A[i*N+j]*x[j]);
            }
            if(residual<fabs(y[i]-x[i])) residual=fabs(y[i]-x[i]);
        }
    };
    void set(double _a){
        G_S::set(_a);
    };
    double* iteration(){
        bool key=false;
        short flag=1;
        while(residual>=Up_Limit_Of_Error && flag>0){
            flag++;
            residual=0;
            if(key){
                for(int i=0;i<N;i++){
                    y[i]=B[i];
                    for(int j=0;j<N;j++)
                        if(i!=j) y[i]+=(A[i*N+j]*x[j]);
                    if(residual<fabs(y[i]-x[i])) residual=fabs(y[i]-x[i]);
                }
            }
            else{
                for(int i=0;i<N;i++){
                    x[i]=B[i];
                    for(int j=0;j<N;j++)
                        if(i!=j) x[i]+=(A[i*N+j]*y[j]);
                    if(residual<fabs(y[i]-x[i])) residual=fabs(y[i]-x[i]);
                }
            }
            key=!key;
        }
        std::cout<<"flag="<<flag<<std::endl;
        return (key?x:y);
    };
};

class SOR : public Jacobi {
    private:
    double omega;
    public:
    SOR(unsigned _n) : Jacobi(_n) {omega=1.9;};
    void set(double** _a){
        Jacobi::set(_a);
    };
    void set(double _a){
        Jacobi::set(_a);
    };
    void set(double* _a){
        for(int i=0;i<N;i++){
            x[i]=B[i]=_a[i]/A[i*N+i];
        }
        residual=0;
        for(int i=0;i<N;i++){
            y[i]=B[i];
            for(int j=0;j<N;j++){
                if(i!=j) y[i]+=(A[i*N+j]*x[j]);
            }
        }
        for(int i=0;i<N;i++){
            y[i]+=((omega-1)*(y[i]-x[i]));
            if(residual<fabs(y[i]-x[i])) residual=fabs(y[i]-x[i]);
        }
    };
    double* iteration(){
        unsigned short flag=1;
        bool key=false;
        while(residual>=Up_Limit_Of_Error && flag>0){
            flag++;
            residual=0;
            if(key){
                for(int i=0;i<N;i++){
                    y[i]=B[i];
                    for(int j=0;j<N;j++){
                        if(i<j) y[i]+=(A[i*N+j]*x[j]);
                        else if(i>j) y[i]+=(A[i*N+j]*y[j]);
                    }
                }
                for(int i=0;i<N;i++){
                    y[i]+=((omega-1)*(y[i]-x[i]));
                    if(residual<fabs(y[i]-x[i])) residual=fabs(y[i]-x[i]);
                }
            }
            else{
                for(int i=0;i<N;i++){
                    x[i]=B[i];
                    for(int j=0;j<N;j++){
                        if(i<j) x[i]+=(A[i*N+j]*y[j]);
                        else if(i>j) x[i]+=(A[i*N+j]*x[j]);
                    }
                }
                for(int i=0;i<N;i++){
                    x[i]+=((omega-1)*(x[i]-y[i]));
                    if(residual<fabs(y[i]-x[i])) residual=fabs(y[i]-x[i]);
                }
            }
            key=!key;
        }
        std::cout<<"flag="<<flag<<std::endl;
        return (key?x:y);
    };
};

#else
//DO NOTHING!
#endif